0.1 Model of Th1-Th2 polarization

Here is the model I am considering for the dynamics of T-cells (where \(T_i\) is the dynamics of any T-cell subpopulation; e.g., \(T_1\) are Th1 cells):

\[\begin{equation} \frac{dT_i}{dt} = b_i + \frac{s_i T_i^2}{S_i^2 + T_i^2}\frac{I_{ij}}{I_{ij}+T_j} - m T_i. \end{equation}\]

The parameter \(b_i\) captures the baseline T-cell production rate that is independent of the T-cell population itself (it has units of T-cells\(\times\)time\(^{-1}\)), for example the T-cell production that is fuelled by cytokine production by macrophages. The fraction \(s_i T_i^2/(S_i^2+T_i^2)\) models the self-promotion of T-cell production fuelled by a T-cell’s own production of cytokines. The parameter \(s_i\) is the maximum rate of self-stimulated production (it also has units of T-cells\(\times\)time\(^{-1}\)). The parameter \(S_i\) is the half-saturation constant for self-promotion. I use a Hill function with an exponent of two based on the work of both Andy Yates and Ed Schrom. The fraction \(I_ij/(I_ij + T_j)\) is the inhibition of proliferation by other T-cell populations, taking values between 0 and 1. The parameter \(I_ij\) gives the abundance of \(T_j\) cells where \(T_i\) proliferation is at half its maximum (it has units of T-cells). The parameter \(m\) gives the baseline rate of immune cell loss, e.g., due to senescence.

The model is really straightforward, but it’s hard to know what parameter values to choose. To help clarify what actually matters, it is useful to nondimensionalize the model. Define new dimensionless state variables \(t_i = T_i/S_i\) and \(\tau = t\times m\), so T-cell density is measured relative to the density at which proliferation is at half its maximum rate. Define the new dimensionless parameters \(\beta_i = b_i/(m\times S_i)\) (the baseline proliferation rate relative to the loss rate when the T cell population is at the density where self-simulated proliferation is half its maximum rate) and \(\sigma_i = s_i/(m\times S_i)\) (the maximum self-stimulated proliferation rate relative to the loss rate at the density where self-simulated proliferation is half its maximum rate). While these parameters are not, in an of themselves, particularly enlightening as to parameter values, their ratio is actually informative: \(\beta_i/\sigma_i\) gives the ratio of baseline to maximum self-stimulated T-cell production. Thus it is reasonable to assume that \(\sigma_i >> \beta_i\). Define the dimensionless parameter \(\iota_{ij} = I_ij/S_i\) (the ratio of the T-cell abundance where cross-inhibition it at half its maximum to the T-cell abudnance where self-stimulation is at half its maximum). Again the value of \(\iota_{ij}\) is informative: a value greater than one implies that cross-inhibition is stronger than self-stimulation, whereas a value less than one implies the opposite.

This gives the dimensionless model

\[\begin{equation} \frac{dt_i}{d\tau} = \beta_i + \frac{\sigma_i t_i^2}{1 + t_i^2}\frac{\iota_{ij}}{\iota_{ij} + t_j} - t_i. \end{equation}\]

The deterministic dynamics of this model can be visualized by plotting the null clines and the deterministic flows. Equilibria occur at the intersections of the clines, and their stability can be visualized from the flow arrow, but I will label each point with a solid triangle if it is a locally stable node, an open triangle if it is an unstable node, and an open diamond if it is an unstable saddle.

Fig. 1 shows the results of this analysis for the case where \((\sigma_1 = \sigma_2 = 2, \beta_1=\beta_2=0.12, \iota_{12}=\iota_{21} = 10)\), implying that self-stimulation is about 20-fold faster than background stimulation, and cross-inhibition is ten-fold stronger than self-stimulation. From this you can see that there are four stable equilibria, corresponding to a resting state (no activation), where both Th1 and Th2 densities are low; a stable co-stimulation state, where both Th1 and Th2 densities are high; and two polarized states where one density is high and the other is low. Thus, the system is multistable and the dynamics will depend on the initial state of the system (deterministically), but also on the initial stochastic dynamics of the system.

## -------------------------------------------------------------------------------
## phaseR: Phase plane analysis of one- and two-dimensional autonomous ODE systems
## -------------------------------------------------------------------------------
## 
## v.2.1: For an overview of the package's functionality enter: ?phaseR
## 
## For news on the latest updates enter: news(package = "phaseR")
Null clines and equilibrium for the dimensionless model.

Figure 1: Null clines and equilibrium for the dimensionless model.

To see what I mean, it is necessary to transform back to the natural scale so that Th1 and Th2 cell abundances are measured at a reasonable scale (right now, they are measured relative to the abundance needed for self-stimulation). A stochastic simulation of the model above would lead to stochastic extinction. We just need to keep our dimensional parameter choices so that they agree with the dimensionless parameters. We will assume that \(m = 1\), to make things simple. Then \(b_i = \beta_i \times S_i\) and \(s_i = \sigma_i \times S_i\), while \(I_{ij} = \iota_{ij} \times S_i\). Thus \(S_i\) (the abundance where self-promotion is half its maximum) becomes the key parameter that defines everything else. Just to make some progress, I will define this as 1000 for both \(T_1\) and \(T_2\).

Fig. 2 shows an example where changing the initial condition of the system leads to really different outcomes for the final state of the system. Starting at (Th1=600, Th2=500), the system approaches a polarized Th1 response (blue lines). Starting at (Th1=500, Th2=600), the system approaches a polarized Th2 response (aqua lines). Starting at (Th1=Th2=550), the system approaches the co-stimulation equilibrium (purple lines).

Null clines, equilibria, and stochastic trajectories from three different initial conditions.

Figure 2: Null clines, equilibria, and stochastic trajectories from three different initial conditions.

However, with the stochastic model, you can also get different outcomes from the same starting point. Here we start in balance, but the different runs give different outcomes (Fig. 3), from Th1 polarized (purple) to Th2 polarized (blue), to co-stimulation (aqua). Exploration of the dynamics of the Th1 and Th2 timecourse suggests that the early response tips the system towards one attractor or the other (Fig. 4).

This indicates that, in addition to understanding how changing the identity of the host or parasite genotype changes the configuration of the isoclines, it is also essential to know where each individual infection “starts”, and to track the initial dynamics of the immune response.

Null clines, equilibria, and three qualitatively different stochastic trajectories with identical initial conditions.

Figure 3: Null clines, equilibria, and three qualitatively different stochastic trajectories with identical initial conditions.

Dynamics of the Th1 and Th2 immune responses from the simulations above.

Figure 4: Dynamics of the Th1 and Th2 immune responses from the simulations above.

0.2 Immune-parasite model

Now we want to add parasites to the model above. If we think about what’s going on biologically, with parasite immune stimulation and immunomodulation, stimulation is mediated through antigen presentation and the cytokine environment, and modulation is mediated through the cytokine environment. Basically, the parasite is able to manipulate the immune response by producing enough cytokine mimics that the cytokine environment tells the host to produce a Th1 response, even though the antigens being presented should trigger a Th2 response. It seems to me that this should be a separate term in the immune dynamic equations, e.g.,: \[\begin{equation} \frac{dT_i}{dt} = b_i + c_i(P) + \frac{s_i T_i^2}{S_i^2 + T_i^2}\frac{I_{ij}}{I_{ij}+T_j} - m T_i. \end{equation}\] where \(c_i(P)\) is the stimulation of Th-\(i\) cells by the parasite.

To get some sense of how this will affect the dynamics of the immune response, notice that if we ignore the actual dynamics of the parasite, we can simply look at how the isoclines change as the sum \(b_i + c_i(P)\) increases, since \(b_i\) will be constant and \(c_i(P)\) will increase (or decrease) as the biomass of parasites increases (or decreases). If both the Th1 and Th2 proliferation rates are increased by the same amount, the resulting configuration of the isoclines is given in Fig. 5 (compare to Fig. 1). The low-activation equilibrium is gone (which makes sense, given that the parasite is there, causing simulation of both Th1 and Th2 responses), but the two polarized responses and the co-stimulation response persist. It is also clear that as the parasite continues to increase in abundance, as long as it continues stimulating both responses equally, the system will settle on a high co-stimulation response. (Note that this also points to a problem with the set of parameters in Fig. 1: the co-stimulation equilibrium is always stable, meaning that even if the parasite is driven extinct, the system will be stuck at this high co-expression equilibrium.) We can deal with this by assuming a much smaller value for baseline immune proliferation, which leads to a situation where the only stable equilibrium is the low activation equilibrium.

derivs = function(t, y, params) {
  sig1 <- params["sig1"]
  sig2 <- params["sig2"]
  beta1 <- params["beta1"]
  beta2 <- params["beta2"]
  i12 <- params["i12"]
  i21 <- params["i21"]
  
  t1 <- y[1]
  t2 <- y[2]
  
  dt1 <- beta1 + sig1*t1^2/(1+t1^2) * i12/(i12+t2) - t1
  dt2 <- beta2 + sig2*t2^2/(1+t2^2) * i21/(i21+t1) - t2
  
  list(c(dt1, dt2))
}

params = c(sig1=2, sig2=2, beta1=0.15, beta2=0.15, i12=10, i21=10)

flows <- flowField(derivs, 
                   xlim=c(0,2),
                   ylim=c(0,2),
                   parameters=params,
                   system="two.dim",
                   state.names=c("Th1","Th2"),
                   add=FALSE)
clines <- nullclines(derivs, 
                     xlim=c(0,2),
                     ylim=c(0,2),
                     parameters=params,
                     col=c(1,2),
                     lwd=2,
                     state.names=c("Th1","Th2"))

eq1 <- findEquilibrium(derivs, y0=c(1.2,1.2), parameters=params, system="two.dim", plot.it=TRUE, summary=FALSE)
eq2 <- findEquilibrium(derivs, y0=c(0.6,1.2), parameters=params, system="two.dim", plot.it=TRUE, summary=FALSE)
eq3 <- findEquilibrium(derivs, y0=c(1.2,0.6), parameters=params, system="two.dim", plot.it=TRUE, summary=FALSE)
eq5 <- findEquilibrium(derivs, y0=c(0.1,0.4), parameters=params, system="two.dim", plot.it=TRUE, summary=FALSE)
eq6 <- findEquilibrium(derivs, y0=c(0.4,0.1), parameters=params, system="two.dim", plot.it=TRUE, summary=FALSE)
The isocline arrangement when parasites are included as a fixed effect on the production rate of both Th1 and Th2 cytokines.

Figure 5: The isocline arrangement when parasites are included as a fixed effect on the production rate of both Th1 and Th2 cytokines.

0.2.1 Stochastic simulations

What I’m going to do now is a stochastic simulation where the parasite abundance grows logistically, with no effect on it by the immune response, just to observe the dynamics of the Th1 and Th2 response. The dynamics of two immune populations are given by \[\begin{equation} \frac{dT_i}{dt} = b1 + c1 P + \frac{s_i T_i^2}{S_i^2 + T_i^2}\frac{I_{ij}}{I_{ij}+T_j} - m T_i. \end{equation}\] As the parasite population grows, the configuration of the \(T_1\) and \(T_2\) isoclines changes. By plotting both the changing isocline configuration and the Th1/Th2 trajectory, I can get a lovely visual of how the parasite population growth warps the flow field through the phase plane, and how that warping affects the dynamics of Th1/Th2 response. In this first simulation, I will assume that the parasite stimulates both immune responses equally well, and that in the absence of infection (\(P=0\)), the only stable equilibrium is the low co-activation equilibrium. The movie shows the dynamics.

derivs = function(t, y, params) {
  s1 <- params["s1"]
  s2 <- params["s2"]
  b1 <- params["b1"]
  b2 <- params["b2"]
  I12 <- params["I12"]
  I21 <- params["I21"]
  S1 <- params["S1"]
  S2 <- params["S2"]
  m <- params["m"]
  c1 <- params["c1"]
  c2 <- params["c2"]
  P <- params["P"]
  
  T1 <- y[1]
  T2 <- y[2]
  
  dT1 <- b1 + c1*P + s1*T1^2/(S1^2+T1^2) * I12/(I12+T2) - m*T1
  dT2 <- b2 + c2*P + s2*T2^2/(S2^2+T2^2) * I21/(I21+T1) - m*T2
  
  list(c(dT1, dT2))
}

gillespie_sim <- function(tmax, y0, params, seed=NULL) {
  if (!is.null(seed))
    set.seed(seed)
  s1 <- params["s1"]
  s2 <- params["s2"]
  b1 <- params["b1"]
  b2 <- params["b2"]
  I12 <- params["I12"]
  I21 <- params["I21"]
  S1 <- params["S1"]
  S2 <- params["S2"]
  m <- params["m"]
  c1 <- params["c1"]
  c2 <- params["c2"]
  bp <- params["bp"]
  Kp <- params["Kp"]
  
  T1 <- y0[1]
  T2 <- y0[2]
  P <- y0[3]
  
  t <- 0
  
  ## set up storage for everything
  out <- array(0, dim=c(1e5, 4))
  colnames(out) <- c("Time", "Th1", "Th2", "P")
  out[1,] <- c(t, T1, T2, P)
  i <- 2
 
  while (t < tmax) {
    ## compute event rates
    prod1 <- b1 + c1*P + s1*T1^2/(S1^2+T1^2) * I12/(I12+T2) 
    prod2 <- b2 + c2*P + s2*T2^2/(S2^2+T2^2) * I21/(I21+T1)
    death1 <- m*T1
    death2 <- m*T2
    birthP <- bp*P*(1-P/Kp)
    rates <- c(prod1, prod2, death1, death2, birthP)
    
    ## what time does the event happen?
    dt <- rexp(1, rate=sum(rates))
    
    ## update t 
    t <- t + dt
    
    ## "wheel of fortune"
    wheel <- cumsum(rates)/sum(rates)
    
    ## which event happens? Draw a random uniform to determine
    rand <- runif(1)
    
    ## if event==1, a new Th1 cell is produced
    ## if event==2, a new Th2 cell is produced
    ## if event==3, a Th1 cell is destroyed
    ## if event==4, a Th2 cell is destroyed
    ## if event==5, a parasite is "born"
    event <- 1 + sum(rand > wheel)
    if (event==1)
      T1 <- T1 + 1
    else if (event==2)
      T2 <- T2 + 1
    else if (event==3)
      T1 <- T1 - 1
    else if (event==4)
      T2 <- T2 - 1
    else 
      P <- P+1
    
    out[i,] <- c(t, T1, T2, P)
    i <- i + 1
    if (i > nrow(out)) ## add more rows
      out <- rbind(out, array(0, dim=c(1e5, 4)))
  }
  out <- out[1:(i-1),]
  return(out)
}


library(magrittr)
set.seed(1234)
seeds <- runif(50, 1, 1000) %>% floor

params = c(S1=1000, S2=1000, s1=2000, s2=2000, b1=0.1, b2=0.1, I12=10000, I21=10000, m=1, c1=1, c2=1, bp=0.2, Kp=200)
out <- gillespie_sim(60, y0=c(5, 5, 10), params=params, seed=seeds[1])

## Plot the isoclines and the Th1/Th2 trajectory as the parasite population grows
## What I want to do is create a series of plots (which will be turned into movies) that shows the Th1/Th2 balance and the isocline configuration every time P crosses a threshold (e.g., every 1 new parasite)
inds <- sapply(seq(10,199,1), function(p) min(which(out[,"P"] > p)))

library(animation)
if(!file.exists("movie2.mp4"))
  saveVideo({ 
    ani.options(interval=0.1)
    for (i in inds) {
      #print(i)
      params2 <- c(params, out[i,"P"])
      flows <- flowField(derivs, 
                         xlim=c(0,2000),
                         ylim=c(0,2000),
                         parameters=params2,
                         system="two.dim",
                         state.names=c("Th1","Th2"),
                         add=FALSE)
      clines <- nullclines(derivs, 
                           xlim=c(0,2000),
                           ylim=c(0,2000),
                           parameters=params2,
                           col=c(1,2),
                           lwd=2,
                           state.names=c("Th1","Th2"))
      lines(out[1:i,2:3],lwd=2, col="orange")
    }}, video.name="movie2.mp4", other.opts = "-pix_fmt yuv420p -b 300k")

You can see that the parasite growth moves faster than the immune dynamics, so the isocline configuration changes faster than the immune response can keep up with it. Because of the fact that the immune dynamics are always relatively slow, there is no opportunity for anything to happen except that the immune dynamics chase a stable co-expression equilibrium as it moves. Presumably, if the parasite growth was a bit slower (or if the immune response reacted a bit faster), there might be an opportunity for the immune dynamics to cross into the basin of attraction for a polarized response. In this first model, I will try slowing the parasite growth down by reducing the value of \(b_p\) (from 0.2 to 0.02).

library(magrittr)
set.seed(1234)
seeds <- runif(50, 1, 1000) %>% floor

params = c(S1=1000, S2=1000, s1=2000, s2=2000, b1=0.1, b2=0.1, I12=10000, I21=10000, m=1, c1=1, c2=1, bp=0.02, Kp=200)
out <- gillespie_sim(500, y0=c(5, 5, 10), params=params, seed=seeds[1])

## Plot the isoclines and the Th1/Th2 trajectory as the parasite population grows
## What I want to do is create a series of plots (which will be turned into movies) that shows the Th1/Th2 balance and the isocline configuration every time P crosses a threshold (e.g., every 1 new parasite)
inds <- sapply(seq(10,199,1), function(p) min(which(out[,"P"] > p)))

library(animation)
if(!file.exists("movie3.mp4"))
  saveVideo({ 
    ani.options(interval=0.1)
    for (i in inds) {
      #print(i)
      params2 <- c(params, out[i,"P"])
      flows <- flowField(derivs, 
                         xlim=c(0,2000),
                         ylim=c(0,2000),
                         parameters=params2,
                         system="two.dim",
                         state.names=c("Th1","Th2"),
                         add=FALSE)
      clines <- nullclines(derivs, 
                           xlim=c(0,2000),
                           ylim=c(0,2000),
                           parameters=params2,
                           col=c(1,2),
                           lwd=2,
                           state.names=c("Th1","Th2"))
      lines(out[1:i,2:3],lwd=2, col="orange")
    }}, video.name="movie3.mp4", other.opts = "-pix_fmt yuv420p -b 300k")

```